Statistical analysis

Litter decomposition

Author

Marcelo Araya-Salas & Andrea Vincent

Published

March 1, 2023

Source code and data found at https://github.com/maRce10/litter_decomposition_EFFEX

Purpose

  • Evaluate role of nutrient availability on litter decomposition

 

Code
cols <- viridis(10, alpha = 0.7)
fill_color <- viridis(10)[7]

# brms models
chains <- 4
iters <- 40000
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"))

# set ggplot2 them
ggplot2::theme_set(theme_classic(base_size = 20))


# standard error
se <- function(x) sd(x)/sqrt(length(x))

1 Nutrient change

1.1 Nitrogen

Content

Code
nutr <- read.csv("./data/raw/litter-nutrients-mas.csv")

nutr$plot.f <- as.factor(nutr$plot)

nutr$days.sc <- scale(nutr$days)

nutr$litter.n.content.prop.initial <- nutr$litter.n.content.perc.initial/100

# remove plot 9
sub.nutr <- nutr[nutr$plot != 9, ]

agg_n <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, mean)

agg_n$sd <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, sd)$litter.n.content.perc.initial

agg_n$se <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, se)$litter.n.content.perc.initial

agg_n$treat <- factor(agg_n$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_n, aes(x = days, y = litter.n.content.perc.initial, color = treat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = litter.n.content.perc.initial +
    se, ymin = litter.n.content.perc.initial - se), width = 0, position = pd) +
    geom_line(size = 1.2, position = pd) + scale_color_viridis_d(alpha = 0.5) +
    labs(x = "Time (days)", y = "Litter N content (% initial)", color = "Treatment") +
    scale_x_continuous(breaks = unique(agg_n$days), labels = unique(agg_n$days)) +
    theme(legend.position = c(0.9, 0.8))

Code
ggsave("./output/Litter_N_content.tiff", dpi = 300)

Download image

Code
fit.n <- brm(litter.n.content.prop.initial ~ treat * days.sc + (1 |
    plot.f), data = nutr, chains = chains, cores = chains, family = Beta(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    prior = prior, file = "./data/processed/litter.n.content.rem_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/litter.n.content.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.2 litter.n.content.rem_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) litter.n.content.prop.initial ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 19631.06 22594.55 272268972
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN -0.340 -0.634 -0.048 1 27796.92 26514.82
treatNP -0.051 -0.344 0.238 1 27673.41 24838.46
treatP -0.060 -0.349 0.231 1 28095.99 25183.40
days.sc -0.503 -0.712 -0.307 1 19631.06 22594.55
treatN:days.sc -0.094 -0.369 0.183 1 25338.01 27759.09
treatNP:days.sc 0.039 -0.227 0.314 1 25217.60 27081.45
treatP:days.sc 0.048 -0.220 0.323 1 23691.70 28287.99

Concentration

Code
agg_conc_n <- aggregate(perc.n ~ colecta + treat + days, nutr, mean)

agg_conc_n$sd <- aggregate(perc.n ~ colecta + treat + days, nutr,
    sd)$perc.n

agg_conc_n$se <- aggregate(perc.n ~ colecta + treat + days, nutr,
    se)$perc.n

agg_conc_n$treat <- factor(agg_conc_n$treat, levels = c("C", "N",
    "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_conc_n, aes(x = days, y = perc.n, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc.n + se, ymin = perc.n -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter N concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_conc_n$days),
    labels = unique(agg_conc_n$days)) + theme(legend.position = c(0.3,
    0.8))

Code
ggsave("./output/Litter_N_concentration.tiff", dpi = 300)

Download image

Code
nutr$prop.n <- nutr$perc.n/100

fit.perc.n <- brm(prop.n ~ treat * days.sc + (1 | plot.f), data = nutr,
    chains = chains, core = chains, family = Beta(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), prior = prior,
    file = "./data/processed/litter.n.perc_model", file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/litter.n.perc_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.3 litter.n.perc_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.n ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 16686.78 23193.79 757288462
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN -0.097 -0.240 0.046 1 25471.79 25380.38
treatNP -0.052 -0.192 0.087 1 25656.85 26419.36
treatP -0.025 -0.163 0.114 1 24461.38 24258.37
days.sc 0.134 0.041 0.227 1 16686.78 23193.79
treatN:days.sc -0.075 -0.205 0.052 1 20777.70 25856.35
treatNP:days.sc -0.030 -0.160 0.099 1 20665.11 26156.06
treatP:days.sc -0.012 -0.138 0.116 1 20899.42 26866.06

1.4 Phosphorus

Content

Code
nutr$litter.p.content.prop.initial <- nutr$litter.p.content.perc.initial/100

# excluding plot 9
agg_p <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, mean)

agg_p$sd <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, sd)$litter.p.content.perc.initial

agg_p$se <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, se)$litter.p.content.perc.initial


agg_p$treat <- factor(agg_p$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_p, aes(x = days, y = litter.p.content.perc.initial, color = treat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = litter.p.content.perc.initial +
    se, ymin = litter.p.content.perc.initial - se), width = 0, position = pd) +
    geom_line(size = 1.2, position = pd) + scale_color_viridis_d(alpha = 0.5) +
    labs(x = "Time (days)", y = "Litter P content (% initial)", color = "Treatment") +
    scale_x_continuous(breaks = unique(agg_p$days), labels = unique(agg_p$days)) +
    theme(legend.position = c(0.9, 0.8))

Code
ggsave("./output/Litter_P_content.tiff", dpi = 300)

Download image

Code
fit.p <- brm(litter.p.content.prop.initial ~ treat * days.sc + (1 |
    plot.f), data = nutr, chains = chains, prior = prior, cores = chains,
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/litter.p.content.rem_model", file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/litter.p.content.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.5 litter.p.content.rem_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) litter.p.content.prop.initial ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 13618.61 17352.52 2063397055
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN -0.090 -1.258 1.081 1 13662.76 18828.78
treatNP 0.009 -1.148 1.172 1 13618.61 17352.52
treatP 0.754 -0.399 1.907 1 14549.30 18770.67
days.sc -0.134 -0.473 0.205 1 17312.08 23553.03
treatN:days.sc 0.001 -0.459 0.459 1 20983.83 27133.63
treatNP:days.sc 0.021 -0.441 0.482 1 21015.75 25715.30
treatP:days.sc -0.168 -0.620 0.288 1 20970.54 26482.58

Concentration

Code
agg_conc_p <- aggregate(perc.p ~ colecta + treat + days, sub.nutr,
    mean)

agg_conc_p$sd <- aggregate(perc.p ~ colecta + treat + days, sub.nutr,
    sd)$perc.p

agg_conc_p$se <- aggregate(perc.p ~ colecta + treat + days, sub.nutr,
    se)$perc.p

agg_conc_p$treat <- factor(agg_conc_p$treat, levels = c("C", "N",
    "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_conc_p, aes(x = days, y = perc.p, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc.p + se, ymin = perc.p -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter P concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_conc_p$days),
    labels = unique(agg_conc_p$days)) + theme(legend.position = c(0.9,
    0.8))

Code
ggsave("./output/Litter_P_concentration.tiff", dpi = 300)

Download image

Code
# convert to proportions to use beta distribution
sub.nutr$prop.p <- sub.nutr$perc.p/100

fit.perc.p <- brm(prop.p ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, cores = chains, family = Beta(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), prior = prior,
    file = "./data/processed/litter.p.perc_model", file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/litter.p.perc_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.6 litter.p.perc_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.p ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 17993.25 23361.17 1797840914
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN -0.106 -0.634 0.409 1 28569.54 28049.17
treatNP 0.031 -0.475 0.530 1 27422.98 28348.78
treatP 0.077 -0.468 0.610 1 28087.50 28163.62
days.sc -0.034 -0.388 0.302 1 17993.25 23361.17
treatN:days.sc -0.001 -0.486 0.478 1 22840.20 28837.19
treatNP:days.sc 0.040 -0.430 0.503 1 21964.87 27252.32
treatP:days.sc 0.028 -0.475 0.531 1 23659.98 28583.96

Takeaways

  • Litter P content is significantly lower in plus N treatment plots than in control plot, after accounting for variation explained by time

 

1.7 Remaining litter

Download image

Code
fit <- brm(prop.litter.rem ~ trat * days.sc + (1 | plot.f), data = dat,
    chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model",
    file_refit = "on_change")

# dat$days.fc <- as.numeric(as.factor(dat$days)) # monotonic
# effect of time fit_mo <- brm(prop.litter.rem ~ trat *
# mo(days.fc) + (1 | plot.f), data = dat, chains = chains,
# family = Beta(), iter = iters, control =
# list(adapt_delta=0.99, max_treedepth=15), cores = chains,
# prior = prior, file =
# './data/processed/prop.litter.rem_model_monotonic', file_refit
# = 'on_change')
Code
extended_summary(read.file = "./data/processed/prop.litter.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.8 prop.litter.rem_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.litter.rem ~ trat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 17979.85 23910.9 1631599382
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tratN -0.038 -0.203 0.125 1 30540.12 29362.46
tratNP 0.074 -0.085 0.233 1 29883.05 29157.97
tratP 0.032 -0.127 0.192 1 30065.88 30156.88
days.sc -0.683 -0.800 -0.570 1 17979.85 23910.90
tratN:days.sc -0.047 -0.207 0.113 1 23277.78 30456.01
tratNP:days.sc 0.106 -0.045 0.259 1 22364.84 28848.92
tratP:days.sc 0.056 -0.097 0.212 1 22771.51 27130.32

1.9 Remaining wood

Download image

 

Code
dat$prop.wood.rem <- dat$perc.wood.rem/100

fit2 <- brm(prop.wood.rem ~ trat * days.sc + (1 | plot.f), data = dat,
    chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.wood.rem_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/prop.wood.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.10 prop.wood.rem_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.wood.rem ~ trat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 16872.71 23764.8 1828689869
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tratN 0.118 -0.128 0.366 1 25420.93 24370.70
tratNP 0.254 0.012 0.496 1 24798.43 26473.85
tratP 0.197 -0.045 0.438 1 25725.64 25244.62
days.sc -0.910 -1.094 -0.737 1 16872.71 23764.80
tratN:days.sc -0.006 -0.244 0.234 1 20661.18 27810.64
tratNP:days.sc 0.097 -0.140 0.333 1 20685.17 28019.73
tratP:days.sc 0.175 -0.058 0.412 1 20811.62 27760.94

2 K

2.1 Litter

Correlation between Silvia’s and Andrea’s K

Code
k_vals <- read.csv("./data/raw/k-values-corr.csv")

cor(k_vals$av.k.litter, k_vals$sil.k.litter)
[1] 0.9327804
Code
cor(k_vals$av.k.wood, k_vals$sil.k.wood)
[1] 0.9621473

Comparing all treatments vs control

Code
agg_kvals <- aggregate(sil.k.litter ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ Treatment, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ Treatment, k_vals, sd)$sil.k.litter
agg_kvals$n <- aggregate(sil.k.litter ~ Treatment, k_vals, length)$sil.k.litter
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 

# composed box plot
ggplot(k_vals, aes(x = Treatment, y = sil.k.litter)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the right
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
   labs(y = "Decomposition constant (k)") +
  # ylim(c(-0.39, 0.145)) +
  geom_text(data = agg_kvals, aes(y = rep(0.4, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
    scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
                              "NP" = "+NP", "P" = "+P"))

Code
ggsave("./output/litter_decomposition_k_by_treatment.tiff", dpi = 300)

Download image

Code
fit_k.litter <- brm(sil.k.litter ~ Treatment + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_model",
    file_refit = "on_change")
Code
extended_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "",
    highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_litter_model.rds")

2.2 k_litter_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.litter ~ Treatment + (1 | quadrat) 20000 4 1 10000 5 0 22486.82 24586.15 224004760
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
TreatmentN 0.238 -0.176 0.659 1 22486.82 24586.15
TreatmentNP -0.156 -0.572 0.261 1 25034.00 26992.34
TreatmentP -0.006 -0.425 0.412 1 23833.14 25611.23

Phosphorus vs no-phosphorus

Code
fit_k.litter.p.np <- brm(sil.k.litter ~ p.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, , control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_p_nop_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/k_litter_p_nop_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.3 k_litter_p_nop_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.litter ~ p.treat + (1 | quadrat) 20000 4 1 10000 8 0 38432.94 27572.89 1236972411
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
p.treatwith.p -0.201 -0.492 0.092 1 38432.94 27572.89

Code
agg_kvals <- aggregate(sil.k.litter ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ p.treat, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ p.treat, k_vals, sd)$sil.k.litter
agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(sil.k.litter ~ p.treat, k_vals, length)$sil.k.litter
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 
k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))

# composed box plot
rn_p_litter <- ggplot(k_vals, aes(x = p.treat, y = sil.k.litter)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_p_litter

Code
ggsave("./output/litter_phosphorus_decomposition_k.tiff", dpi = 300)

Nitrogen vs no-nitrogen

Code
fit_k.litter.n.nn <- brm(sil.k.litter ~ n.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_n_no_n_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/k_litter_n_no_n_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.4 k_litter_n_no_n_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.litter ~ n.treat + (1 | quadrat) 20000 4 1 10000 17 0 37023.21 27505.83 956380015
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
n.treatwith.n 0.045 -0.254 0.345 1 37023.21 27505.83

Code
agg_kvals <- aggregate(sil.k.litter ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ n.treat, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ n.treat, k_vals, sd)$sil.k.litter
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(sil.k.litter ~ n.treat, k_vals, length)$sil.k.litter
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))

# composed box plot
rn_n_litter <- ggplot(k_vals, aes(x = n.treat, y = sil.k.litter)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_n_litter

Code
ggsave("./output/litter_nitrogen_decomposition_k.tiff", dpi = 300)

2.5 Wood

Comparing all treatments vs control

Code
agg_kvals <- aggregate(sil.k.wood ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ Treatment, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ Treatment, k_vals, sd)$sil.k.wood
agg_kvals$n <- aggregate(sil.k.wood ~ Treatment, k_vals, length)$sil.k.wood
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 

# composed box plot
ggplot(k_vals, aes(x = Treatment, y = sil.k.wood)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
   labs(y = "Decomposition constant (k)") +
  # ylim(c(-0.39, 0.145)) +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
    scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
                              "NP" = "+NP", "P" = "+P"))

Code
ggsave("./output/wood_decomposition_k_by_treatment.tiff", dpi = 300)

Download image

Code
fit_k.wood.p.treat <- brm(sil.k.wood ~ Treatment + (1 | quadrat),
    data = k_vals, chains = chains, family = gaussian(), iter = iters *
        2, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/k_wood_treatment_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/k_wood_treatment_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.6 k_wood_treatment_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.wood ~ Treatment + (1 | quadrat) 80000 4 1 40000 23 0 96744.72 104808.8 528819192
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
TreatmentN 0.034 -0.537 0.610 1 96900.55 107769.7
TreatmentNP -0.490 -1.064 0.087 1 96744.72 104808.8
TreatmentP -0.469 -1.044 0.105 1 97055.26 107363.6

3 Phosphorus vs no-phosphorus

Code
agg_kvals <- aggregate(sil.k.wood ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ p.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ p.treat, k_vals, sd)$sil.k.wood


agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(sil.k.wood ~ p.treat, k_vals, length)$sil.k.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 


k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))


# composed box plot
rn_p_wood <- ggplot(k_vals, aes(x = p.treat, y = sil.k.wood)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_p_wood

Code
ggsave("./output/phosphorus_decomposition_k.tiff", dpi = 300)

Download image

Code
fit_k.wood.p.no.p <- brm(sil.k.wood ~ p.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_wood_p_no_p_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/k_wood_p_no_p_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

3.1 k_wood_p_no_p_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.wood ~ p.treat + (1 | quadrat) 20000 4 1 10000 4 0 31952.25 26238.54 326631910
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
p.treatwith.p -0.496 -0.892 -0.1 1 31952.25 26238.54

4 Nitrogen vs no-nitrogen

Code
agg_kvals <- aggregate(sil.k.wood ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ n.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ n.treat, k_vals, sd)$sil.k.wood
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(sil.k.wood ~ n.treat, k_vals, length)$sil.k.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))

# composed box plot
rn_n_wood <- ggplot(k_vals, aes(x = n.treat, y = sil.k.wood)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_n_wood

Code
ggsave("./output/nitrogen_decomposition_k.tiff", dpi = 300)

Download image

Code
fit_k.wood.n.no.n <- brm(sil.k.wood ~ n.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_wood_n_no_n_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/k_wood_n_no_n_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

4.1 k_wood_n_no_n_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.wood ~ n.treat + (1 | quadrat) 20000 4 1 10000 7 0 30908.88 27135.72 1683655054
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
n.treatN 0.008 -0.414 0.426 1 30908.88 27135.72

5 Nutrient content by remaining litter mass

5.1 Nitrogen

Code
nutr$prom.hoja.reman.sc <- scale(nutr$prom.hoja.reman)

ggplot(data = nutr, aes(x = prom.hoja.reman, y = perc.n, color = treat)) +
    geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter N concentration (%)",
    color = "Treatment") + scale_color_viridis_d(alpha = 0.5) + geom_smooth(method = "lm",
    se = FALSE) + scale_x_reverse()

Code
ggsave("./output/nitrogen_by_litter_mass.tiff", dpi = 300)

Download image

Code
fit_pern_by_rem_leave <- brm(perc.n ~ prom.hoja.reman.sc * treat +
    (1 | plot), data = nutr, chains = chains, family = gaussian(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/n_per_by_leave_rem_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/n_per_by_leave_rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

5.2 n_per_by_leave_rem_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) perc.n ~ prom.hoja.reman.sc * treat + (1 | plot) 20000 4 1 10000 0 0 22712.6 26086.89 202030287
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
prom.hoja.reman.sc -0.152 -0.265 -0.038 1 22712.60 27352.64
treatN -0.161 -0.345 0.025 1 29404.33 26086.89
treatNP -0.044 -0.224 0.138 1 29226.14 26624.22
treatP -0.037 -0.218 0.145 1 28414.57 26316.85
prom.hoja.reman.sc:treatN 0.060 -0.095 0.212 1 27004.70 29358.59
prom.hoja.reman.sc:treatNP -0.103 -0.269 0.064 1 28381.31 29425.85
prom.hoja.reman.sc:treatP -0.013 -0.174 0.150 1 27640.20 30237.33

5.3 Phosphorus

Code
ggplot(data = sub.nutr, aes(x = prom.hoja.reman, y = perc.p, color = treat)) +
    geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter P concentration (%)",
    color = "Treatment") + scale_color_viridis_d(alpha = 0.5) + geom_smooth(method = "lm",
    se = FALSE) + scale_x_reverse()

Code
ggsave("./output/phosphorus_by_litter_mass.tiff", dpi = 300)

Download image

Code
sub.nutr$prom.hoja.reman.sc <- scale(sub.nutr$prom.hoja.reman)

fit_perp_by_rem_leave <- brm(perc.p ~ prom.hoja.reman.sc * treat +
    (1 | plot), data = sub.nutr, chains = chains, family = gaussian(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/p_per_by_leave_rem_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/p_per_by_leave_rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

5.4 p_per_by_leave_rem_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) perc.p ~ prom.hoja.reman.sc * treat + (1 | plot) 20000 4 1 10000 0 0 18469.62 19550.38 8818990
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
prom.hoja.reman.sc 0.001 -0.002 0.004 1 21082.94 25940.57
treatN -0.005 -0.012 0.002 1 18785.39 19550.38
treatNP 0.000 -0.007 0.007 1 18469.62 20718.11
treatP 0.003 -0.004 0.011 1 18490.36 20587.53
prom.hoja.reman.sc:treatN 0.000 -0.004 0.005 1 24965.09 28102.51
prom.hoja.reman.sc:treatNP -0.002 -0.006 0.003 1 25491.25 27830.25
prom.hoja.reman.sc:treatP -0.001 -0.005 0.004 1 26931.48 27838.27

Takeaways

  • Nitrogen percentage, but no Phosphorus percentage, increases along with remaining litter mass

 

6 Nutrient ratios

6.1 C:N

Code
# excluding plot 9
agg_cn <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, mean)

agg_cn$sd <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, sd)$cn.mol.kg

agg_cn$se <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, se)$cn.mol.kg


agg_cn$treat <- factor(agg_cn$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_cn, aes(x = days, y = cn.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = cn.mol.kg + se, ymin = cn.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter C:N ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_cn$days),
    labels = unique(agg_cn$days)) + theme(legend.position = c(0.9,
    0.8))

Code
ggsave("./output/cn_ratio_through_time.tiff", dpi = 300)

Download image

Code
fit.cn <- brm(cn.mol.kg ~ treat * days.sc + (1 | plot.f), data = nutr,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/cn_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/cn_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.2 cn_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) cn.mol.kg ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 23616.14 27188.36 206225990
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN 0.309 -1.373 2.015 1 33345.57 27188.36
treatNP -0.189 -1.881 1.509 1 31802.40 27824.88
treatP -0.182 -1.863 1.514 1 30821.46 28785.15
days.sc -3.101 -4.204 -1.996 1 23616.14 27323.41
treatN:days.sc 0.751 -0.752 2.245 1 29469.69 30884.24
treatNP:days.sc 0.439 -1.062 1.942 1 29505.97 31955.60
treatP:days.sc 0.111 -1.391 1.614 1 29297.05 31633.29

6.3 N:P

Code
# excluding plot 9
agg_np <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    mean)

agg_np$sd <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    sd)$np.mol.kg

agg_np$se <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    se)$np.mol.kg


agg_np$treat <- factor(agg_np$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_np, aes(x = days, y = np.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = np.mol.kg + se, ymin = np.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter N:P ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_np$days),
    labels = unique(agg_np$days)) + theme(legend.position = c(0.2,
    0.8))

Code
ggsave("./output/np_ratio_through_time.tiff", dpi = 300)

Download image

Code
fit.np <- brm(np.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/np_model",
    file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/np_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.4 np_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 5.3) np.mol.kg ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 23138.62 23571.24 226626333
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN 1.211 -3.076 5.572 1 27662.79 23728.89
treatNP -2.401 -6.601 1.896 1 26490.49 23571.24
treatP -3.135 -7.696 1.500 1 27059.33 23906.00
days.sc 4.407 1.755 7.072 1 23138.62 26390.49
treatN:days.sc -1.765 -5.391 1.885 1 26580.63 27306.41
treatNP:days.sc -1.780 -5.403 1.901 1 27640.74 28593.42
treatP:days.sc -1.415 -5.318 2.486 1 28596.79 29590.75

6.5 C:P

Code
# excluding plot 9
agg_cp <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    mean)

agg_cp$sd <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    sd)$cp.mol.kg

agg_cp$se <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    se)$cp.mol.kg

agg_cp$treat <- factor(agg_cp$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_cp, aes(x = days, y = cp.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = cp.mol.kg + se, ymin = cp.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter C:P ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_cp$days),
    labels = unique(agg_cp$days)) + theme(legend.position = c(0.9,
    0.8))

Code
ggsave("./output/cp_ratio_through_time.tiff", dpi = 300)

Download image

Code
fit.cp <- brm(cp.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, family = gaussian(), iter = iters, cores = chains,
    prior = prior, file = "./data/processed/cp_model", file_refit = "on_change")
Code
extended_summary(read.file = "./data/processed/cp_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.6 cp_model

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 106.3) cp.mol.kg ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 30235.23 26541.32 1179254390
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatN 0.386 -19.206 20.050 1 30235.23 27580.58
treatNP -0.159 -19.492 19.226 1 31461.97 26541.32
treatP -0.300 -19.898 19.294 1 32550.30 27705.93
days.sc -6.177 -21.408 9.060 1 30959.55 27830.18
treatN:days.sc 0.319 -17.923 18.305 1 32235.63 27315.59
treatNP:days.sc -4.800 -23.038 13.518 1 32491.62 27863.26
treatP:days.sc -2.648 -20.973 15.782 1 32066.82 27837.57


7 Combined plots

7.1 Remaining mass

Code
agg_rem_w$substrate <- "Wood"
agg_rem$substrate <- "Litter"

names(agg_rem) <- names(agg_rem_w) <- c("colecta", "trat", "days",
    "perc.rem", "sd", "se", "substrate")

agg_rem_pooled <- rbind(agg_rem, agg_rem_w)

ggplot(agg_rem_pooled, aes(x = days, y = perc.rem, color = trat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = perc.rem +
    se, ymin = perc.rem - se), width = 0, position = pd) + geom_line(size = 1.2,
    position = pd) + scale_color_viridis_d(alpha = 0.5, labels = c("Control",
    "+N", "+P", "+NP")) + labs(x = "Time (days)", y = "Remaining mass (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_rem_pooled$days),
    labels = unique(agg_rem_pooled$days)) + facet_wrap(~substrate,
    nrow = 2)

Code
ggsave("./output/remaining_mass_combined.tiff", dpi = 300)

Download image

7.2 Litter content

Code
agg_p$nutrient <- "P"
agg_n$nutrient <- "N"

names(agg_n) <- names(agg_p) <- c("colecta", "treat", "days", "perc",
    "sd", "se", "nutrient")

agg_np2 <- rbind(agg_n, agg_p)

ggplot(agg_np2, aes(x = days, y = perc, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc + se, ymin = perc -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N",
        "+P", "+NP")) + labs(x = "Time (days)", y = "Litter nutrient content (% initial)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_p$days),
    labels = unique(agg_p$days)) + facet_wrap(~nutrient, nrow = 2)

Code
ggsave("./output/litter_content_combined.tiff", dpi = 300)

Download image

7.3 Concentration

Code
agg_conc_p$nutrient <- "P"
agg_conc_n$nutrient <- "N"

names(agg_conc_n) <- names(agg_conc_p) <- c("colecta", "treat", "days",
    "perc", "sd", "se", "nutrient")

agg_conc_np <- rbind(agg_conc_n, agg_conc_p)

ggplot(agg_conc_np, aes(x = days, y = perc, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc + se, ymin = perc -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N",
        "+P", "+NP")) + labs(x = "Time (days)", y = "Litter nutrient concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_conc_p$days),
    labels = unique(agg_conc_p$days)) + facet_wrap(~nutrient, nrow = 2,
    scales = "free_y")

Code
ggsave("./output/concentration_combined.tiff", dpi = 300)

Download image

7.4 Nutrient ratios

Code
agg_cn$nutrients <- "C:N"
agg_cp$nutrients <- "C:P"
agg_np$nutrients <- "N:P"

names(agg_cn) <- names(agg_cp) <- names(agg_np) <- c("colecta", "treat",
    "days", "mol.kg", "sd", "se", "nutrients")

agg_ratios <- rbind(agg_cn[, names(agg_cn)], agg_cp[, names(agg_cn)],
    agg_np[, names(agg_cn)])

ggplot(agg_ratios, aes(x = days, y = mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = mol.kg + se, ymin = mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N",
        "+P", "+NP")) + labs(x = "Time (days)", y = "Litter nutrient ratios",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_ratios$days),
    labels = unique(agg_ratios$days)) + facet_wrap(~nutrients, nrow = 3,
    scales = "free_y")

Code
ggsave("./output/nutrient_ratios_combined.tiff", dpi = 300, width = 10,
    height = 12)

Download image

7.5 Decomposition constant

Code
rn_n_wood_dat <- rn_n_wood$data[, c("n.treat", "sil.k.wood")]
rn_n_litter_dat <- rn_n_litter$data[, c("n.treat", "sil.k.litter")]
rn_p_wood_dat <- rn_p_wood$data[, c("p.treat", "sil.k.wood")]
rn_p_litter_dat <- rn_p_litter$data[, c("p.treat", "sil.k.litter")]

names(rn_p_wood_dat) <- names(rn_n_wood_dat)  <- names(rn_p_litter_dat) <- names(rn_n_litter_dat) <- c("treatment", "k")
rn_p_wood_dat$substrate <- rn_n_wood_dat$substrate <- "wood"
rn_p_litter_dat$substrate <- rn_n_litter_dat$substrate <- "litter"

rn_n_wood_dat$nutrient <- rn_n_litter_dat$nutrient <- "N"
rn_p_litter_dat$nutrient <- rn_p_wood_dat$nutrient <- "P"

rn_k_dat <- rbind(rn_n_wood_dat, rn_p_wood_dat, rn_n_litter_dat, rn_p_litter_dat)                             
                             
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  # geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
    facet_grid(~ substrate)

Code
ggsave("./output/decomposition_k_combined_v1.tiff", dpi = 300, width = 10, height = 12)

Download image

Code
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  # geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
    facet_grid(substrate ~ nutrient, scales = "free_x")

Code
ggsave("./output/decomposition_k_combined_v12.tiff", dpi = 300, width = 10, height = 12)

Download image

8 Combined model diagnostics

Code
check_rds_fits(path = "./data/processed", html = TRUE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) cn.mol.kg ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 11186.547 18465.237 206225990
2 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 106.3) cp.mol.kg ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 3800.096 3310.208 1179254390
3 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.litter ~ Treatment + (1 | quadrat) 20000 4 1 10000 5 0 10254.769 11806.319 224004760
4 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.litter ~ n.treat + (1 | quadrat) 20000 4 1 10000 17 0 8997.572 8327.074 956380015
5 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.litter ~ p.treat + (1 | quadrat) 20000 4 1 10000 8 0 10806.876 12607.636 1236972411
6 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.wood ~ n.treat + (1 | quadrat) 20000 4 1 10000 7 0 8165.879 9671.253 1683655054
7 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.wood ~ p.treat + (1 | quadrat) 20000 4 1 10000 4 0 8190.692 11336.204 326631910
8 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) sil.k.wood ~ Treatment + (1 | quadrat) 80000 4 1 40000 23 0 37124.811 54565.753 528819192
9 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) litter.n.content.prop.initial ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 12453.012 17014.307 272268972
10 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.n ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 14603.282 16880.616 757288462
11 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) litter.p.content.prop.initial ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 8457.580 15360.879 2063397055
12 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.p ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 16204.903 17372.163 1797840914
13 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) perc.n ~ prom.hoja.reman.sc * treat + (1 | plot) 20000 4 1 10000 0 0 12014.609 18513.789 202030287
14 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 5.3) np.mol.kg ~ treat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 11291.391 17616.689 226626333
15 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5) perc.p ~ prom.hoja.reman.sc * treat + (1 | plot) 20000 4 1 10000 0 0 7081.653 10836.794 8818990
16 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.litter.rem ~ trat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 13040.090 18549.848 1631599382
17 b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) prop.wood.rem ~ trat * days.sc + (1 | plot.f) 20000 4 1 10000 0 0 11933.976 17029.142 1828689869

 

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggdist_3.2.0      brms_2.18.0       Rcpp_1.0.9        viridis_0.6.2    
 [5] viridisLite_0.4.1 ggridges_0.5.4    posterior_1.3.1   ggrepel_0.9.1    
 [9] cowplot_1.1.1     tidybayes_3.0.1   ggplot2_3.4.0     readxl_1.3.1     
[13] knitr_1.42        kableExtra_1.3.4 

loaded via a namespace (and not attached):
  [1] backports_1.4.1      systemfonts_1.0.4    plyr_1.8.7          
  [4] igraph_1.3.5         splines_4.1.0        svUnit_1.0.6        
  [7] crosstalk_1.2.0      TH.data_1.1-0        rstantools_2.2.0    
 [10] inline_0.3.19        digest_0.6.31        htmltools_0.5.4     
 [13] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
 [16] remotes_2.4.2        RcppParallel_5.1.5   matrixStats_0.62.0  
 [19] xts_0.12.2           sandwich_3.0-1       svglite_2.1.0       
 [22] prettyunits_1.1.1    colorspace_2.0-3     rvest_1.0.3         
 [25] textshaping_0.3.5    xfun_0.36            dplyr_1.0.10        
 [28] callr_3.7.3          crayon_1.5.2         jsonlite_1.8.4      
 [31] lme4_1.1-27.1        survival_3.2-11      zoo_1.8-11          
 [34] glue_1.6.2           gtable_0.3.1         emmeans_1.8.1-1     
 [37] webshot_0.5.4        distributional_0.3.1 pkgbuild_1.4.0      
 [40] rstan_2.21.7         abind_1.4-5          scales_1.2.1        
 [43] mvtnorm_1.1-3        DBI_1.1.1            xaringanExtra_0.7.0 
 [46] miniUI_0.1.1.1       xtable_1.8-4         stats4_4.1.0        
 [49] StanHeaders_2.21.0-7 DT_0.26              htmlwidgets_1.5.4   
 [52] httr_1.4.4           threejs_0.3.3        arrayhelpers_1.1-0  
 [55] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1.9000      
 [58] farver_2.1.1         utf8_1.2.2           labeling_0.4.2      
 [61] tidyselect_1.2.0     rlang_1.0.6          reshape2_1.4.4      
 [64] later_1.3.0          munsell_0.5.0        cellranger_1.1.0    
 [67] tools_4.1.0          cli_3.6.0            generics_0.1.3      
 [70] evaluate_0.20        stringr_1.5.0        fastmap_1.1.0       
 [73] ragg_1.1.3           yaml_2.3.7           processx_3.8.0      
 [76] purrr_1.0.0          packrat_0.9.0        pbapply_1.6-0       
 [79] nlme_3.1-152         mime_0.12            projpred_2.0.2      
 [82] formatR_1.11         xml2_1.3.3           compiler_4.1.0      
 [85] bayesplot_1.9.0      shinythemes_1.2.0    rstudioapi_0.14     
 [88] gamm4_0.2-6          sketchy_1.0.2        tibble_3.1.8        
 [91] stringi_1.7.12       ps_1.7.2             Brobdingnag_1.2-9   
 [94] lattice_0.20-44      Matrix_1.5-1         nloptr_1.2.2.2      
 [97] markdown_1.3         shinyjs_2.1.0        tensorA_0.36.2      
[100] vctrs_0.5.2          pillar_1.8.1         lifecycle_1.0.3     
[103] bridgesampling_1.1-2 estimability_1.4.1   httpuv_1.6.6        
[106] R6_2.5.1             promises_1.2.0.1     gridExtra_2.3       
[109] codetools_0.2-18     boot_1.3-28          colourpicker_1.2.0  
[112] MASS_7.3-54          gtools_3.9.3         assertthat_0.2.1    
[115] rprojroot_2.0.3      withr_2.5.0          shinystan_2.6.0     
[118] multcomp_1.4-17      mgcv_1.8-36          parallel_4.1.0      
[121] gghalves_0.1.3       grid_4.1.0           minqa_1.2.4         
[124] tidyr_1.1.3          coda_0.19-4          rmarkdown_2.20      
[127] shiny_1.7.3          base64enc_0.1-3      dygraphs_1.1.1.6